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ABSTRACT 

We discuss a model for treating chemical enrichment by Type II and Type la super- 
nova (SNII and SNIa) explosions in simulations of cosmological structure formation. 
Our model includes metal-dependent radiative cooling and star formation in dense 
collapsed gas clumps. Metals are returned into the diffuse interstellar medium by star 

■ particles using a local smoothed particle hydrodynamics (SPH) smoothing kernel. A 

variety of chemical abundance patterns in enriched gas arise in our treatment owing 
to the different yields and lifetimes of SNII and SNIa progenitor stars. In the case of 
SNII chemical production, we adopt metal-dependent yields. Because of the sensitive 
dependence of cooling rates on metallicity, enrichment of galactic haloes with metals 

■ can in principle significantly alter subsequent gas infall and the build up of the stellar 
| components. Indeed, in simulations of isolated galaxies we find that a consistent treat- 
, ment of metal-dependent cooling produces 25 per cent more stars outside the central 

region than simulations with a primordial cooling function. In the highly-enriched cen- 
tral regions, the evolution of baryons is however not affected by metal cooling, because 
here the gas is always dense enough to cool. A similar situation is found in cosmolog- 
ical simulations because we include no strong feedback processes which could spread 
metals over large distances and mix them into unenriched diffuse gas. We demonstrate 
this explicitly with test simulations which adopt suprasolar cooling functions leading 
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to large changes both in the stellar mass and in the metal distributions. We also find 
that the impact of metallicity on the star formation histories of galaxies may depend 
on their particular evolutionary history. Our results hence emphasise the importance 
of feedback processes for interpreting the cosmic metal enrichment. 



Key words: methods: numerical - galaxies: abundances - galaxies: formation - galax- 
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1 INTRODUCTION 

Over the last decades, our knowledge of the chemical prop- 
erties of the Universe and, in particular, of galaxies, has 
improved dramatically. Observations of the Local Universe 
(e.g. Garnett & Shields 1987; Skillman, Kennicutt & Hodge 
1989; Brodie & Huchra 1991; Zaritsky, Kennicutt & Huchra 
1994; Mushotzky et al. 1996; Ettori et al. 2002; Tremonti et 
al. 2004; Lamareille et al. 2004) as well as at intermediate 
and high redshifts (e.g. Prochaska & Wolfe 2002; Adelberger 
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et al. 2003; Kobulnicky et al. 2003; Lilly, Carollo & Stockton 
2003; Shapley et al. 2004) have resulted in a quite detailed 
picture of the chemical history of the stellar populations and 
of the interstellar and intergalactic media. It is an important 
challenge for the theory of galaxy formation to explain these 
observational results in the context of the current cosmolog- 
ical paradigm. Numerical simulations are an important tool 
to make specific predictions for galaxy formation theories 
and to confront them with observations, because they can 
accurately track the growth of structure in the dark matter 
and baryonic components. It is hence natural to ask what 
such simulations predict for the chemical properties of the 
universe. 

Currently, hydrodynamic cosmological simulations in- 
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eluding radiative cooling and star formation, at best, 
coarsely reproduce the observed properties of galaxies. They 
fail when scrutinised in detail. Among the most prominent 
problems are "catastrophic angular momentum loss" (e.g. 
Navarro & Benz 1991; Navarro & White 1994) and the dif- 
ficulty in finding a consistent modelling of supernova (SN) 
feedback (e.g. Navarro & White 1993; Metzler & Evrard 
1994; Yepes et al. 1997; Marri & White 2003; Springel & 
Hernquist 2003, hereafter SH03). The first issue mainly af- 
fects galaxy morphology, size and kinematics, while the lat- 
ter primarily influences the efficiency of star formation. Both 
are the subject of intense research. Reproducing the chemi- 
cal properties of observed galaxies in detail may be viewed 
as an intertwined third problem. 

Regarding the angular momentum problem, recent sim- 
ulations have shown some promising advances. Dommguez- 
Tenreiro, Tissera & Saiz (1998) found that the formation of 
compact stellar bulges can stabilize the discs and thereby 
prevent the angular momentum loss during violent minor 
merger events, such that disc systems that resemble observed 
galaxies are obtained. Abadi et al. (2003) have also demon- 
strated the importance of a dense, slowly rotating spheroidal 
component, and pointed out its relevance for comparing sim- 
ulated and observed galaxies consistently. More recently, 
Robertson et al. (2004) were able to produce large, sta- 
ble disc systems in their cosmological simulations which are 
comparable in size to spiral galaxies. In their approach, they 
introduced a subresolution model for star formation and 
feedback which pressurizes the star-forming gas and sta- 
bilizes discs against fragmentation. We note that all these 
studies agree on the need for a self-consistent treatment of 
SN feedback as a crucial mechanism to regulate star forma- 
tion and to reproduce discs similar to observational coun- 
terparts. 

SN explosions are thought to play a fundamental role in 
the evolution of galaxies since they are considered the most 
efficient and ubiquitous mechanism for the ejection of metals 
and energy into the interstellar medium (ISM). Both chem- 
ical and energy feedback can affect the condensation of gas 
and consequently the evolution of galactic systems. On one 
hand, the presence of metals in the ISM affects gas cooling 
times because the radiative cooling rate depends sensitively 
on metallicity (Sutherland & Dopita 1993, hereafter SD93). 
On the other hand, the release of energy is crucial for reg- 
ulating star formation through the heating and disruption 
of cold gas clouds, and for producing outflows which can 
transport enriched material into the intergalactic medium 
(e.g. Lehnert & Heckman 1996; Dahlem, Weaver & Heckman 
1998; Frye, Broadhurst & Bem'tez 2002; Rupke, Veilleux & 
Sanders 2002; Martin 2004). As suggested by observations 
and theory, the largest outflows (e.g. Larson 1974; White & 
Rees 1978; Dekel & Silk 1986; White & Frenk 1991) should 
be able to develop in small systems because of their shal- 
lower potential wells. Because in hierarchical galaxy forma- 
tion scenarios large systems are formed by the aggregation 
of smaller systems, energy feedback from SNe is then ex- 
pected to have important effects for nearly all systems in 
the different stages of galaxy formation. 

Modelling of chemical feedback has been addressed al- 
ready in numerous studies, mostly with the aim to repro- 
duce the chemical properties of certain types of galaxies, or 
of clusters of galaxies (e.g. Larson 1976; Tinsley & Larson 



1979; Burkert & Hensler 1988; White & Frenk 1991; Burkert, 
Truran & Hensler 1992; Ferrini et al. 1992; Theis, Burkert 
& Hensler 1992; Steinmetz & Miiller 1994; Kauffmann 1996; 
Chiappini, Matteucci & Gratton 1997; Kauffmann & Char- 
lot 1998; Boisser & Prantzos 2000; Valdarnini 2003). As first 
showed by Mosconi et al. (2001), in order to properly model 
the process of metal enrichment in galaxy formation it is nec- 
essary to consider the full cosmological growth of structure 
in which mergers and interactions have important effects on 
the star formation process and the dynamical evolution of 
galaxies (see also Kawata & Gibson 2003; Tornatore et al. 
2004; Okamoto et al. 2005). 

Perhaps the most fundamental motivation for the rele- 
vance of metal enrichment for galaxy formation is based on 
the fact that the cooling rate of baryons depends on metal- 
licity (SD93). Radiative cooling in turn allows the conden- 
sation of gas into dense and cold clouds, which form the 
reservoir of material available for the formation of stars. As 
chemical enrichment is a result of star formation, a chem- 
ical feedback cycle emerges, in which metals can in prin- 
ciple significantly accelerate the transformation of baryons 
into stars. In hierarchical scenarios these processes depend 
significantly on the chemical prescription adopted (Kaellan- 
der & Hultman 1998; Kay et al. 2000). These affect galaxy 
properties such as the luminosity function (White & Frenk 
1991). For these reasons, a detailed analysis of the effects 
of metal-dependent cooling is clearly important for galaxy 
formation studies. Substantial numerical challenges in sim- 
ulations of galaxy formation come from the large dynamic 
range required, and from the multiphase character of the 
ISM. Gas in very different physical states coexists in the 
ISM, but the standard formulation of smoothed particle hy- 
drodynamics (SPH, e.g. Gingold & Monaghan 1977; Lucy 
1977) is not well suited to deal with this situation. A number 
of attempts have been made to solve this problem, ranging 
from an ad-hoc decoupling of phases (e.g. Pearce et al. 1999, 
2001; Marri & White 2003) to simple analytic sub-resolution 
models (SH03) that work with an effective equation of state. 

Our new model for SN feedback in SPH cosmological 
simulations has been divided into two stages: chemical and 
energy feedback. In this work we address the description 
of the chemical feedback, which is based on the previous 
model of Mosconi et al. (2001). In a forthcoming paper, we 
will present the second part of this work: the implementa- 
tion of energy feedback in the framework of a multiphase 
scheme for the ISM, which extends the approach of Marri 
& White (2003). While a treatment of chemical enrichment 
without energy feedback is an incomplete picture, we think 
the complexity of the problem merits a step- wise discussion 
of our model such that the primary physical effects can be 
better understood. 

This paper is organized as follows. In Section^ we sum- 
marize the numerical implementation of the chemical model, 
and in Section |H] we analyse the dependence of the results 
on numerical resolution and input parameters, using sim- 
ulations of isolated galaxies. We then study in Section 2] 
the properties of galaxies formed in cosmological simula- 
tions and examine how they depend on the chemical model. 
Finally, in Section Qj] we give our conclusions. 



© 0000 RAS, MNRAS 000, [I]-?? 



Chemical feedback in SPH 3 



2 NUMERICAL IMPLEMENTATION 

Our new model for the production and ejection of chem- 
ical elements by SNe is based in part on the approach in 
Mosconi et al. (2001). We have implemented our scheme 
in the TreePM/SPH code GADGET-2, an improved version 
of the public code GADGET (Springel, Yoshida & White 
2001) which manifestly conserves energy and entropy where 
appropriate (Springel & Hernquist 2002). Note that we do 
not use the multiphase treatment and the feedback model 
developed by SH03 for this code, but we do include their 
treatment of ultraviolet (UV) background. Also, we use a 
different parameterisation for star formation, as described 
in detail below. 

Chemical elements are synthesized in stellar interiors 
and ejected into the ISM by SN explosions. In order to model 
the production and distribution of metals in the context of 
the simulations, there are three ingredients we have to con- 
sider: the SN rate (i.e. number of exploding stars per time 
unit), the chemical yields (i.e. chemical material ejected in 
explosions) and the typical lifetimes of stars, which deter- 
mine the characteristic time of the metal release. 

In our model, we include a separate treatment of Type 
II and Type la SNe (SNII and SNIa). These two types of 
SNe originate from different stellar populations, and have 
different rates, yields and typical time-scales. We use dif- 
ferent yields for SNII and SNIa which, in the case of SNII, 
are metal-dependent. Hence, we will treat SNII and SNIa 
separately, owing to the need to model their different char- 
acteristics. For example, SNII produce most of the chemi- 
cal elements, except for iron, which is mainly produced by 
SNIa. Another difference between these two types of SN is 
that SNII are the endpoint of the evolution of massive stars 
with short lifetimes, in contrast to SNIa, which result from 
the evolution of binary systems with lifetimes of ~ 1 Gyr. 

We assume that initially gas particles have primordial 
abundances: Xh = 0.76 and Yho = 0.24, and we consider 
the enrichment by the following elements, assuming produc- 
tion only by SNII and SNIa: H, 4 He, 12 C, 16 , 24 Mg, 28 Si, 
56 Fe, 14 N, 20 Ne, 32 S, 40 Ca and 62 Zn. Note that we do not 
consider nucleosynthesis by intermediate-mass stars. Conse- 
quently, we will concentrate our analysis on elements such as 
O or Fe where restricting to SN production may be a good 
approximation. In the rest of this Section, we describe the 
most important aspects of the chemical model, namely the 
star formation and cooling prescriptions, and our scheme for 
the production and ejection of metals. 



2.1 Star formation 

We assume that gas particles are eligible for star forma- 
tion if they are denser than a critical value (p > = 
7.0 x 10 _26 g cm -3 , where p denotes gas density) and lie 
in a convergent flow (divw < 0) 1 . For these particles, we 
assume a star formation rate (SFR) per unit volume equal 
to 
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where c is a star formation efficiency (we adopt c = 0.1) and 
Tdyn = l/^/ArirGp is the dynamical time of the particle. We 
base the creation of new stellar particles on the stochastic 
approach of SH03 (see also Lia, Portinari & Carraro 2002). 
To this end, each gas particle eligible for star formation in 
a given time-step is assigned a probability p* of forming a 
star particle given by 

i c At\ 



p* = 



expl 



(2) 



where At is the integration time-step of the code, m is the 
mass of the gas particle, and m* is defined as m* = mo/N g 
with mo being the original mass of gas particles at the begin- 
ning of the simulation and N g a parameter which determines 
the number of 'generations' of stars formed from a given gas 
particle (we assume N g = 2). If a random number drawn 
from a uniform distribution in the unit interval is smaller 
than p«, we form a new stellar particle of mass m» and re- 
duce the mass of the gas particle accordingly. Once the gas 
particle mass has become smaller than mo/N g , we instead 
turn it into a star particle. Newly born stars are assumed 
to have the same element abundances as the gas mass from 
which they form. 

The main advantage of this stochastic approach is that 
all particles are either purely stellar or purely gas, which 
avoids the use of 'hybrid' particles with both gaseous and 
stellar components. The latter would artificially force the 
stellar component to evolve dissipately like the gas. Note 
that, as a consequence of the star formation scheme, the 
number of baryonic particles in the simulation does not re- 
main constant. Unlike SH03, we also exchange the masses 
of heavy elements explicitly between gas and star particles. 
As a result, the mass spectrum of stellar and gas particles is 
slightly washed out around otherwise sharp discrete values. 

2.2 Chemical production and distribution 

Our numerical implementation of the enrichment model has 
three main components: SNII element production, SNIa el- 
ement production and metal distribution. 

SNII element production 

In order to estimate the SNII rate we use a Salpeter Initial 
Mass Function (IMF) with lower and upper mass cut-offs of 
0.1 and 40 Mq, respectively, and assume that stars more 
massive than 8Mq end their lives as SNII. For the chem- 
ical production, we adopt the yields of Woosley & Weaver 
(1995). These are metal-dependent yields and vary differ- 
ently for different elements. We use half of the iron yield 
of WW95, as it is often adopted (e.g. Timmes, Woosley & 
Weaver 1995). For the sake of simplicity, in most of our sim- 
ulations we have assumed that SNII explode (i.e. produce 
and eject the chemical elements) within an integration time- 
step of the code (which is typically < 10 6 yr), because SNII 
originate in massive stars which have very short character- 
istic lifetimes of the order of 10 7 yr. However, we have also 
analysed the effect of relaxing this instantaneous recycling 
approximation (IRA) for SNII (see Section [3.21 . 

SNIa element production 

As a progenitor model, we adopt the W7 model of Thiele- 
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mann, Nomoto & Hashimoto (1993), which assumes that 
SNIa explosions originate from CO white dwarf systems in 
which mass is transferred from the secondary to the primary 
star until the Chandrasekhar mass is exceeded and an explo- 
sion is triggered. It is generally assumed that the lifetime of 
such a binary system is in the range rsNia =[0.1, 1 ] Gyr, de- 
pending on the age of the secondary star (Greggio 1996) . In 
our model, the material produced by a SNIa event is ejected 
when a time rsNia after the formation of the exploding star 
has elapsed. We choose this time randomly within a given 
range (in most of our experiments we assume rsNia = [0.1, 1] 
Gyr). In order to estimate the number of SNIa, we adopt an 
observationally motivated range for the relative ratio of SNII 
and SNIa rates (e.g. van den Bergh 1991). We use the chemi- 
cal yields of Thieleman et al. (1993) for the metal production 
itself. 

Metal distribution 

We distribute chemical elements ejected in SN explo- 
sions within the gaseous neighbours of exploding star par- 
ticles, using the usual SPH kernel interpolation technique 
(Mosconi et al. 2001). Each neighbouring gas particle re- 
ceives a fraction of the ejected metals according to its kernel 
weight. For this purpose, gas particle neighbours of stellar 
particles are identified whenever metal distribution has to 
take place, and smoothing lengths for these stars are cal- 
culated by setting a desired number of gas neighbours that 
should be enclosed in the smoothing length (we have as- 
sumed the same criterium used for the smoothing lengths of 
gas particles). This is particularly important for the enrich- 
ment by SNIa because of the time delay between the forma- 
tion of the stars and the ejection of metals. Note that after 
a star particle is created, it is only subject to gravitational 
forces, unlike the gas particles, which continue to be sub- 
ject to hydrodynamical forces as well. Consequently, as the 
system evolves, star particles can release metals associated 
with SNIa explosions in a different environment from where 
they were born. This, together with the different lifetimes 
and yields of SNII explosions, can then produce non-trivial 
variations in abundance ratios. 



2.3 Gas cooling 

It is well known that the cooling rate of thin astrophysical 
plasmas depends sensitively on metallicity, in such a way 
that at a given density, the cooling rate is higher for a higher 
metallicity. In Fig.Qwe show the cooling function computed 
by SD93 from primordial (i.e. no heavy elements) to supra- 
solar abundance ([Fe/H]=+0.5). As one can see from this 
figure, the cooling function can show very large variations 
with the metal abundance of the gas, depending on the tem- 
perature range considered. For example, the usual estimates 
of the cooling time (r coo i) for a gas particle at the critical 
density (p*) and at a temperature of T» = 4 x 10 4 K yield 
that Tcooi is 50 times larger for primordial gas than for a 
mildly suprasolar medium ([Fe/H]= 0.5), 18 times larger 
than that of a solar abundance one, and still twice as large 
as the corresponding time for gas with [Fe/H]= —1.5. The 
temperature range of [10 5 — 10 6 ] K shows the largest differ- 
ences, right in the temperature range relevant to the abun- 
dant population of dwarf galaxies. 







1 


i i i i i j 

0.5 ^ 


11111111111 




-21 


- 










- 


/ /-05/T 














B 








or. 

In 

<u 


-22 


- 


AT / / Ai>"\ Vs ^ \ 


\ 
< 






' V //a\ 2 -0 \ \ V \ ^\ >■ 


log 






XJ \-3.o\ 






-23 


- 


Primordial\\ 





4 5 6 7 

log (T/K) 



Figure 1. Cooling rates for gas of different metallicities, as given 
by the [Fe/H] abundance, from suprasolar ([Fc/H]=0.5) to pri- 
mordial. The plot is based on data adapted from SD93. 

Including effects owing to metal line cooling is there- 
fore crucial for realistic modelling of galaxy formation be- 
cause it affects directly the relation between cooling time 
and dynamical time of forming galaxies. We therefore in- 
clude a treatment of radiative cooling consistent with the 
metal content of the gas component produced by the en- 
richment model. We use the tables computed by SD93 for 
the metal-dependent cooling functions, and interpolate from 
them for other metallicities as needed. For particles with iron 
abundance larger than 0.5, we adopt the suprasolar cooling 
rate (i.e. that for [Fe/H]> +0.5). In cosmological simula- 
tions, we also include a redshift-dependent photo-heating 
UV background, with an intensity evolution that follows the 
model of Haardt & Madau (1996). 



3 RESULTS FOR ISOLATED GALAXY 
MODELS 

As first tests of our model, we consider simulations of iso- 
lated haloes set up to form disc galaxies at their centre. This 
allows us to examine the dependence of the results on resolu- 
tion and on different parameters of the chemical enrichment 
model. The initial conditions consist of a static dark mat- 
ter potential corresponding to a Navarro, Frenk & White 
(1996, 1997) profile of concentration c = 20, and a baryonic 
gas phase initially in hydrostatic equilibrium in this poten- 
tial. Our typical system has virial mass M 200 = 10 12 ft -1 M s 
(h = 0.7), 10 per cent of which is in the form of baryons 
(Mbar = 10 n /i -1 Mq). The initial radius of the system is 
T200 = 160/i _1 kpc and the gas has an initial angular mo- 
mentum characterized by a spin parameter A = 0.1. We 
have followed the evolution of this system for 1.5 Td yn , where 
i"dyn = 7 Gyr is the dynamical time at r2oo- 

These idealized initial conditions yield a simple model 
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for disc formation, which is an ideal test bench for the per- 
formance and validity of the code, and the dependence of the 
results on the free parameters. We stress however that these 
models are not meant to provide a realistic scenario for the 
whole galaxy formation process. Results for full cosmological 
simulations will be discussed in Section 4, where we analyse 
the properties of galactic objects formed in A cold dark mat- 
ter (CDM) models with virial masses of M vir « 10 12 ft _1 M , 
comparing them also to the results of the idealized tests dis- 
cussed here. 

3.1 Numerical resolution 

Both the star formation process and the distribution of met- 
als can be affected by numerical resolution. Particularly im- 
portant for our model is the fact that the smoothing length 
of the particles (i.e. the radius which encloses « 32 gaseous 
neighbours) determines the region where the newly released 
chemical elements are injected. Mosconi et al. (2001) found 
that chemical abundances could be oversmoothed in low res- 
olution simulations because metals are distributed over a 
comparatively large volume and, consequently, are more ef- 
ficiently mixed with surrounding gas. On the other hand, 
high numerical resolution could also produce spurious re- 
sults if the metal mixing mechanism becomes too inefficient. 
Provided other physical mixing mechanisms are not operat- 
ing, the simulation results might show an artificially large 
dispersion in gas metallicities. 

In order to test the effects of numerical resolution, we 
carried out simulations of the isolated disc system described 
above using three different initial numbers of gas particles: 
2500 (Rl), 10000 (R2) and 40000 (R3), with mass resolu- 
tions of 4 x 10 7 , 10 7 and 0.25 x 10 7 h~ l M Q for Rl, R2 and 
R3, respectively. We have set the gravitational softening to 
0.4 ft -1 kpc. In these three runs we use metal-dependent 
cooling functions for the gas particles according to the tables 
by SD93. In TableQwe summarize the main characteristics 
of these simulations. 

In Fig. |3 we show the evolution of the SFRs obtained 
for Rl (dotted line), R2 (solid line) and R3 (dashed line). 
Most of the stellar mass (92, 88 and 86 per cent in Rl, 
R2 and R3, respectively) is formed at times t < 3 Gyr (43 
per cent of the dynamical time) in the three runs. Only at 
later times, we find that the lowest resolution test simulation 
(Rl) yields a significantly lower result for the SFR, indicat- 
ing that the number of gas particles has dropped so much by 
the conversion into stars that discrete effects become appre- 
ciable. However, the final stellar mass fractions (i.e. stellar 
mass formed divided by the total baryonic mass) are very 
similar in Rl, R2 and R3. 

This can be also appreciated from Fig.|2] where we show 
the mass- weighted metallicity profiles for Rl (dotted line), 
R2 (solid line) and R3 (dashed line), in the inner region 
and for the stellar (left column) and gaseous (right column) 
components and for three different times, t = 0.13 Gyr (up- 
per panels), t = 2.5 Gyr (middle panels) and t = 7 Gyr 
(lower panels). We can see that the stellar abundance pro- 
files are insensitive to resolution. The mean iron abundance 
differs by less than 0.03 dex between the different simula- 
tions. Recall that the stellar metallicity of newly forming 
stars is given by the local gas metallicity. This is directly 
visible at the beginning of the simulation where the gaseous 
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Figure 2. Evolution of the SFR for our test simulations of iso- 
lated disc galaxies, as a function of different numerical resolution. 
Rl (dotted line) has 2500 particles, R2 (solid line) 10000 particles, 
and R3 (dashed line) 40000 particles. 

and stellar components show a similar level of enrichment 
in star- forming regions (upper panels in Fig. |3J, while at 
later stages of the evolution (middle and lower panels) this 
constraint becomes hidden as a result of the superposition 
of stars of different ages. 

However, as can be seen from the right panels of Fig. 
the mass-weighted metallicity profiles for the gas compo- 
nents show stronger differences with resolution than their 
stellar counterparts. At the beginning of the evolution, the 
profiles of the three test simulations are still very similar, 
although the enriched gas is more spread in the highest- 
resolution test. However, the gas profiles differ more at later 
times when factors such as variations in the consumption 
time of the gas and differences in the spatial distribution of 
stars begin to affect the radial distribution of metals. At the 
end of the simulations (lower right panel) , the mean gaseous 
iron abundance for Rl is significantly lower than those for 
the other two simulations. This decrease in the metallicity 
is produced by the infall of less enriched gas. This gas is not 
becoming dense enough to trigger star formation activity at 
the level of R2 and R3, and the infall tends to dilute the 
metal content in the ISM. Hence, the fact that the gas com- 
ponent in Rl shows a lower iron abundance is driven by its 
poorer gas resolution. Note that as a result of the consump- 
tion of most of the gas, the final number of gas particles in Rl 
is very low (~ 60 in the inner 30 kpc) and, as a consequence, 
the distribution of metals in the gas component is strongly 
affected by numerical noise. On the other hand, the results 
for the tests with 10000 and 40000 particles converge quite 
well both in their stellar and gaseous chemical properties. 
Note that the rms scatter of particles within each radial bin 
is similar for the different tests indicating that, on average, 
the metal mixing is not significantly affected by resolution 
if we consider at least an initial number of particles of 2500 
or more. 

3.2 Dependence on assumed supernova 
char acter ist ics 

In this Section we analyse the dependence of our results on 
important input parameters of the chemical model, namely 
the IRA for SNII, and the rate and lifetimes associated with 
SNIa. For this purpose we run three additional test simula- 
tions. In runs PI and P2, we change the SNIa parameters 
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Figure 3. Metallicity profiles for the stellar (left panels) and gaseous (right panels) components of the idealized disc tests performed 
with different resolution: Rl (dotted line, 2500 particles), R2 (solid line, 10000 particles) and R3 (dashed line, 40000 particles). The 
different rows show different times of the evolution: t = 0.13 Gyr (upper panels), t = 2.5 Gyr (middle panels) and t = 7 Gyr (lower 
panels). The error bars correspond to the rms scatter of [Fe/H] around the mean in each radial bin. 



Table 1. Chemical model parameters for the different tests of the idealized isolated galaxy: initial number of gas particles -/V gas , typical 
lifetimes associated with SNI1 (rgNn) an d SNIa (rgNla) in Gyr, rate of SNIa and cooling functions adopted. 
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compared with the simulation R2, which we take as a fiducial 
reference model, while in run P3 we release the IRA condi- 
tion for SNII. The main characteristics of these simulations 
are listed in Tabled 

In Fig. 2] we compare the evolution of the SFR of the 
standard R2 run (solid line) with PI (dotted line), in which 
the SNIa rate was increased. We find no significant differ- 
ences between the SFRs of these tests. However, we find 
a larger difference in their stellar age-metallicity relations 



(AMR), as seen in Fig.|^| although the results appear consis- 
tent within the standard deviation based on counting statis- 
tics. In this plot, an increase in the SNIa rate translates into 
a shift towards larger values of stellar metallicities, since 
the ISM is more strongly enriched by heavy elements. A 
change of the SNIa rate also impacts other chemical prop- 
erties of the systems. For example, in Fig. |S| we show the 
[O/Fe] abundance as a function of [Fe/H] for stars in R2 
and PI at Td yn . The main difference is found for the high- 
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Figure 4. SFR for the experiments of the idealized disc galaxy 
run with different SN parameters: R2 (solid line, standard test), 
PI (dotted line, higher SNIa rate), P2 (dashed line, IRA for SNIa) 
and P3 (dashed-dotted line, relaxing the IRA for SNII). 




2 4 

t [Gyr] 



Figure 5. AMR for the stellar component in the tests of the ideal- 
ized disc galaxy run with different SN parameters: R2 (solid line, 
standard test), PI (dotted line, higher SNIa rate), P2 (dashed 
line, IRA for SNIa), and P3 (dashed-dotted line, relaxing the IRA 
for SNII). The error bars correspond to the rms scatter of [Fe/H] 
around the mean. 

metallicity material ([Fe/H] > -0.5) which shows in PI a 
lower cv-enhancement compared to that in R2, as expected. 

We have also tested the sensitivity of the results on the 
lifetime of binary systems associated with SNIa. The impor- 
tance of including SNIa has been previously pointed out by 
Greggio & Renzini (1983) and Mosconi et al. (2001), among 
others. Motivated by these previous results, we carried out a 
simulation P2 where we assumed an IRA for SNIa explosions 
(see Table [TJ. The main effect of instantaneously releasing 
the chemical production of SNIa can be appreciated from 
Fig. H]] which shows the mean a-enhancement of the simu- 
lated stellar population. These abundance ratios exhibit a 
behaviour which is in open disagreement with observations 
(Pagel 1997). Note, however, that these observational results 
correspond to solar neighbour stars, while the simulated ones 
include information from the whole stellar component. 

Finally, in order to analyse the effects of the IRA for 
SNII, we tested a more detailed calculation of the lifetime 
of massive stars. To this end, we convolved the IMF with the 
lifetime of massive stars at different intervals of metallicity 
(Raiteri, Villata & Navarro 1996). In this way, we obtained a 
range of mean lifetimes varying with stellar metallicity from 
6.3 x 10 6 yr for Z = 5 x 10 _4 Z© to 1.6 x 10 7 yr for Z = Z Q , 
where Z denotes mean stellar metallicity and Zq is the solar 
metallicity. We run a test of our idealized initial condition 




-4 —3 -7 -1 

[Fe/H] 



Figure 6. [O/Fc] abundance as a function of [Fe/H] for the stellar 
component in experiments of the idealized disc galaxy run with 
different SN parameters: R2 (solid line, standard test), PI (dotted 
line, higher SNIa rate), P2 (dashed line, IRA for SNIa), and P3 
(dashed-dotted line, relaxing the IRA for SNII). The error bars 
correspond to the rms scatter of [O/Fe] around the mean. 




t [Gyr] 



Figure 7. SFR for the idealized disc test run with suprasolar (SC, 
dashed line), primordial (PC, dotted line) and metal-dependent 
(R2, solid line) cooling functions. 



with the same chemical parameters used in R2 but including 
these metallicity-dependent lifetimes for SNII (P3, see Ta- 
bic 0. We found no significant differences between the two 
tests, either in the SFR or in the iron abundance, as can 
be seen from Figs. 2] an d The similar behaviour of these 
two runs is also found in other chemical properties of the 
systems, such as the relation shown in Fig. |S| Recall that a 
drawback of a detailed description of chemical enrichment 
can be its high computational cost, so an important prac- 
tical goal is to keep the model sufficiently simple to allow 
its use in large-scale cosmological simulations. While in this 
context the IRA for SNIa in general cannot reproduce cer- 
tain observed trends, for SNII it appears to be a simple and 
sufficiently accurate simplification, at least for these simple 
tests that do not include energy feedback. 



3.3 Effects of metal-dependent cooling 

Finally, we performed two more runs of our idealized iso- 
lated galaxy in order to highlight effects owing to a metal- 
dependent cooling function. In these fiducial tests, we as- 
sumed a constant metallicity, either primordial (PC) or 
suprasolar (SC, [Fe/H]= +0.5) for the cooling function. The 
parameters for star formation and SN rates were the same 
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Figure 8. Face-on (left-hand panels) and edge-on (right-hand panels) surface metallicity distribution for the stellar component of the 
isolated disc galaxy tests run with primordial (PC, upper row), metal-dependent (R2, middle row) and suprasolar (SC, lower row) cooling 
rate functions. The metallicity scale is also shown. 



as in the standard simulation R2 discussed earlier (see Ta- 
ble [TJ. 

In Fig. Q we compare the evolution of the SFR for R2 
(solid line), PC (dotted line) and SC (dashed line). In the 
early stages (t < 4 Gyr), the test run with suprasolar abun- 



dance cooling function (SC) results in a SFR higher by up 
to a factor of 3 compared with the simulation R2 with self- 
consistent cooling function, while the run PC with primor- 
dial cooling lies lower by up to a factor 1.5. After the first 
4 Gyr, the SFR in the SC run decreases significantly, indi- 



© 0000 RAS, MNRAS 000, [T^?? 



Chemical feedback in SPH 9 



eating that most of the gas has already been consumed and 
transformed into stars. Likewise, after this period, simula- 
tion PC has the highest residual level of SFR as a result 
of the larger amount of gas left over and available for star 
formation. 

As a result of the differences in the level of star for- 
mation activity among the runs with different cooling func- 
tions, the metallicity distribution is also affected. In Fig. |H] 
we show metallicity maps for the stellar components corre- 
sponding to face-on (left panels) and edge-on (right panels) 
projections of the discs, for the PC (upper panel), R2 (mid- 
dle panel) and SC (lower panel) simulations. Comparing the 
upper to the lower panels, the distribution of stars is more 
extended and the level of enrichment is higher in the SC 
run. In fact, the increase in the cooling efficiency produces 
an increase in the star formation activity at larger radius, 
increasing the amount of metals and affecting their relative 
distribution. This can also be appreciated from Fig.|5J which 
shows the fraction of baryons in stellar form f(r) as a func- 
tion of radius for tests R2 (solid line), PC (dotted line) and 
SC (dashed line) at Td yn . These fractions were calculated 
by summing up the stellar mass in radial bins, normalized 
to the total baryonic mass within the corresponding bin. 
In particular, at 20 h^ 1 kpc the distribution functions for 
R2 and SC are 25 and 35 per cent higher than that for the 
PC run, respectively. Note that the impact of using different 
cooling functions becomes more important as we go to outer 
regions where, in the limiting case of SC, the high cooling 
efficiency leads to a substantial star formation activity even 
outside 30 h^ 1 kpc. The behaviour found in this figure is 
present already in the early stages of the evolution. The ef- 
ficient transformation of gas into stars at larger radius also 
limits the supply of baryons to inner regions. Hence, the 
whole mass distribution is affected by the use of different 
cooling functions. This is clear from Fig. 1101 where we show 
the density profiles of the stellar (upper panel) and gaseous 
(lower panel) components for these tests. 

Finally, in Fig. [TT] we show the mass-weighted oxygen 
profiles for the stellar (left panels) and gaseous (right pan- 
els) components for the tests PC (dotted line), R2 (solid 
line) and SC (dashed line), at three different times: t = 0.13 
Gyr (upper panels), t = 2.5 Gyr (middle panels), t = 7 Gyr 
(lower panels) . We find no significant differences in the abun- 
dance distributions for the stars in these tests at any time. 
However, the metallicity of the gaseous component is more 
sensitive to the adopted cooling functions. At the beginning 
of the simulations, the results for the three runs are similar 
but after the initial collapse (middle panel) the level of en- 
richment increases from PC to R2 and SC, even though the 
differences remain quite small overall. At Tdyn (lower panel), 
both the PC and R2 simulations show decreasing metallic- 
ity profiles, conversely to the case of SC which shows a flat 
relation. Here metals are more uniformly distributed and 
the ISM is enriched out to larger radii. Recall that the gas 
metallicities reflect an instantaneous state of the chemical 
properties of the system, while the stellar population gives 
an integrated account of the history of the chemical proper- 
ties of the galaxy. 

In summary, the results of this section show that a de- 
tailed study of the metallicity properties of galaxies requires 
a self-consistent treatment of radiative cooling that accounts 
for metal line cooling. The resulting effects will, however, de- 
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Figure 9. Fraction of baryons in stellar form f(r) for the ideal- 
ized disc tests run with suprasolar (SC, dashed line), primordial 
(PC, dotted line) and metal-dependent (R2, solid line) cooling 
functions. 
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Figure 10. Density profiles for the stellar (upper panel) and 
gaseous (lower panel) components for the idealized disc tests run 
with suprasolar (SC, dashed line), primordial (PC, dotted line) 
and metal-dependent (R2, solid line) cooling functions. 

pend strongly on how heavy elements are spread into diffuse 
unenriched gas. 



4 RESULTS FOR COSMOLOGICAL 
SIMULATIONS 

We here discuss first results for full cosmological simula- 
tions of the hierarchical growth of galaxies, allowing us to 
more realistically assess effects of chemical enrichment on 
galaxy properties. To this end, we have run a cosmolog- 
ical simulation (C2) in a CDM universe with cosmologi- 
cal parameters fi m = 0.3, Qa = 0.7, fib — 0.04, Ho = 
100/i _1 km s _1 Mpc _1 , h = 0.7, and cr 8 = 0.9. We used a 
periodic simulation box of 10 ft -1 Mpc on a side, populated 
initially with N = 2 x 80 3 particles, which yields an ini- 
tial mass resolution of « 1.4 x 10 s and ss 2.16 x 10 7 h~ x 
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Figure 11. Oxygen profiles for the stellar (left panel) and gaseous (right panel) components in experiments of the idealized disc galaxy 
run with different cooling functions: R2 (solid line, metal-dependent cooling), PC (dotted line, cooling function for primordial abundance) 
and SC (dashed line, cooling function for suprasolar abundance). The error bars correspond to the rms scatter of 12 +log (O/H) around 
the mean. 



Mq for dark matter and gas particles, respectively. We have 
adopted a maximum gravitational softening of 5 hT 1 kpc 
for dark matter, gas and star particles. We used a metal- 
dependent, self-consistent cooling function for the gas com- 
ponent. Table |21 summarizes the main characteristics of this 
run. 

In order to test the dependence of our results on mixing 
processes and metal-dependent cooling, we have also carried 
out two additional runs with the same parameters, but a 
constant metallicity for the cooling function. As a limiting 
case of a very enriched medium, in simulation C3 we used 
suprasolar abundance (see Table |5J, while in simulation CI 
we took primordial cooling in order to model a situation 
where mixing processes are completely absent. These two 
cases should hence bracket reality. 

In our cosmological simulations, we identified galac- 
tic objects at z = as virialized structures using a 
density-contrast criterion with overdensity 5p/p « 178 Qq ' 6 
(White, Efstathiou & Frenk 1993). In our analysis we focus 
on six haloes with masses similar to the idealized disc test 
R2 (~ 1O 12 M ). They are resolved with a similar or higher 
number of particles (N > 10000, see Table |3J. Note that as 
gas/star particles do not have the same mass, the number 
of particles in these components does not, however, trace 



the mass. The main properties of these objects are given in 
Table H 

In Fig. 1121 we compare the star formation histories for 
the selected galactic objects in CI (dotted line), C2 (solid 
line) and C3 (dashed line), as a function of look-back time, 
i.e. t(z) = 1 — [1 + z]~ 3//2 . Note that the global features of 
the SFRs are preserved in the three runs, although the de- 
tailed levels of activity are different. In general, the SFRs of 
the galactic objects formed in C3 are systematically higher 
than the corresponding ones in C2, as expected. Note that 
because galactic objects formed in C3 have higher SFRs than 
those in C2, their final stellar masses are much larger than 
the corresponding ones in C2 (see Table [IJ. If we compare 
the SFRs obtained for CI and C2, although the differences 
are small, the SFRs in CI are systematically lower than 
those in C2. It is noteworthy that the impact of metallic- 
ity on the star formation history seems to depend on the 
particular evolutionary histories of the systems, as can be 
appreciated from Fig. 1121 

In Fig. 1131 we show the mean stellar mass fraction (f(r)) 
as a function of radius, estimated from the six selected galac- 
tic objects in CI (dotted line), C2 (solid line) and C3 (dashed 
line) . As in Section 13.31 the distribution function for each 
galactic object was calculated by summing up the stellar 
mass in bins of radius normalized to the total baryonic mass 
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Table 2. Chemical model parameters for the cosmological tests: initial number of gas and dark matter particles N gaB and Wdm, typical 
lifetimes associated with SNII (tsnii) and SNIa (rgNla) in Gyr, rate of SNIa and cooling functions adopted. 
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Figure 12. SFR for the six selected galactic objects from the cosmological simulations run with primordial (dotted line), metal-dependent 
(solid line) and suprasolar (dashed lines) cooling functions. 



Table 3. Number of dark matter (JVdm), gas (-/V gas ), and stars (JV s tar) within the virial radius for the six galactic objects (GLOs) 
identified in the cosmological tests run with primordial (CI), metal-dependent (C2) and suprasolar (C3) cooling functions. 
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Table 4. Virial mass (M v ; r ), stellar mass (M star , both in I0 10 h 1 Mq) and optical radius (r opt in h 1 kpc) for the six galactic objects 
(GLOs) identified in the cosmological tests run with primordial (CI), metal-dependent (C2) and suprasolar (C3) cooling functions. 
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within the corresponding bin. As we found for the tests of the 
isolated galaxies, there are clear differences between the re- 
sults for the three runs. Although in the inner 10 hT 1 kpc the 
simulations yield similar stellar fractions, within 40 kpc, 
galactic objects in C2 have transformed a larger fraction of 
the gas into stars compared with those in C3, while galactic 
objects in CI show the opposite behaviour. At larger radii, 
the mean distribution function of stellar mass for the galac- 
tic objects in C3 is significantly larger than its counterpart 
for C2, indicating that as a result of the higher cooling effi- 
ciency, the gas in C3 has been transformed into stars further 
away from the centre of mass of the objects. Conversely, the 
primordial run shows fewer stars in the outer region. Note 
that efficient energy feedback is needed to prevent excessive 
star formation in the outer regions, as it can be observed 
from Fig. [TT] 

In order to assess how these differences affect the chem- 
ical properties of the formed galactic objects, we analysed 
the metal indicator 12 + log(0/H) for the gaseous and stel- 
lar components in CI, C2 and C3. Consistent with our re- 
sults for the idealized disc tests, we found that the metallic- 
ity of the stellar components is not sensitive to the cooling 
function adopted, in contrast to the gaseous components. 
In Fig. 1141 we compare the gaseous oxygen profiles for the 
galaxy samples in CI (dotted lines), C2 (solid lines) and C3 
(dashed lines) . Although there is a trend for systems in C3 to 
have a higher level of enrichment and flatter mass-weighted 
abundance profiles compared to their counterparts in C2 and 
CI, how much metals affect these relations seems to depend 
on the particular history of evolution of each galaxy, as can 
be inferred from Figs. I12l and 1141 

Finally, for a comparison with observational results, we 
performed linear regression fits to the gaseous oxygen pro- 
files of our simulated galactic objects within the optical ra- 
dius. By definition, the optical radius encloses 83 per cent of 
the baryonic stellar mass within the virial radius. In Fig. 1151 
we show the zero-point of these fits as a function of the cor- 
responding slopes for the galactic objects in CI (triangles), 
C2 (diamonds) and C3 (asterisks). We also include obser- 
vational results from Zaritsky et al. (1994, filled circles). 
In general, the properties of the simulated galactic objects 
agree reasonably well with the range obtained in observa- 
tional studies. However, the simulated galactic objects iden- 
tified in C2 tend to have steeper profiles and higher zero- 
points compared with observations. Conversely, galactic ob- 
jects in C3 exhibit flatter profiles. A similar behaviour was 
also found in our tests for the isolated disc galaxies run 
with primordial (PC, squared triangle), metal-dependent 
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Figure 13. Mean stellar mass fraction as a function of radius for 
the six selected objects in the cosmological simulations run with 
primordial (CI, dotted line), metal-dependent (C2, solid line) and 
suprasolar (C3, dashed line) cooling functions as a function of 
radius. 



(R2, squared diamond) and suprasolar (SC, squared aster- 
isk) cooling functions. The fact that these simulations do 
not reproduce the whole observed range of gradients could 
be related to the lack of a treatment of energy feedback 
by SNe. The latter, if realistically modelled, is expected to 
trigger mass outflows from the star-forming regions into the 
halo. Since these regions are highly enriched, metals may be 
transported outwards, affecting the metallicity distributions 
and flattening the metallicity profiles. Our simulation run 
with suprasolar cooling function was performed in order to 
assess the effects of the presence of metals in these outer 
regions on the gas cooling, mimicking an extreme case of 
metal enrichment by energy feedback. 



5 CONCLUSIONS 

We have implemented a new model for chemical evolution 
within the cosmological SPH code GADGET-2. We account 
for nucleosynthesis from SNIa and SNII separately. Radia- 
tive cooling of the gas according to self-consistent metallicity 
values is realized by means of precomputed look-up tables. 
Our implementation was designed to allow a reasonably ac- 
curate treatment of the main features of stellar evolution 
while still not becoming computationally too expensive. 

Using a number of test simulations both of isolated 
galaxies and of cosmological structure formation, we have 
validated our implementation and established some first ba- 
sic results for the effects of chemical enrichment. An impor- 
tant consequence of the presence of metals is the increase in 
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Figure 14. Oxygen profiles for the gaseous components of the six selected galactic objects from the cosmological simulations run with 
primordial (dotted lines), metal-dependent (solid lines) and suprasolar (dashed lines) cooling functions. The error bars correspond to the 
rms scatter of 12 + log(0/H) around the mean. 



the cooling efficiency that it can cause, leading to acceler- 
ated star formation. In our isolated galaxy simulations, the 
SFR for test runs with metal-dependent cooling function 
can be up to a factor of 1.5 higher than the corresponding 
one for a test with a primordial cooling function, and up 
to a factor of 3 lower than the one for suprasolar cooling. 
We found that these differences do not originate in the gas 
in the central regions, because it is already cold and dense 
anyway. Instead, the differences are caused by more diffuse 
gas at outer radii which has yet to cool. Here the increase 
of the cooling efficiency allows condensation of gas and star 
formation to occur further away from the centre of the sys- 
tems. 

An important implication of this result is that strong 
effects owing to metal-dependent cooling can only be ex- 
pected when heavy elements are efficiently transported and 
mixed into diffuse gas that has yet to cool. In our cosmolog- 
ical simulations, which lacked a physical model for energy 
feedback, such a large-scale spreading of metals only occurs 
at a low level. As perhaps was to be expected, we found that 
simulations with metal-dependent cooling produced similar 
results as simulations with a primordial cooling function. We 
expect however that this could be very different if a model 
for efficient metal distribution is added to our simulations. 
An explicit demonstration of this has been provided by our 



test simulation where we assumed suprasolar cooling rate 
functions, resulting in substantially elevated star formation 
activity and different abundance profiles. Note that in this 
paper we explored galaxy scale systems which cool very effi- 
ciently. Larger effects from chemical mixing are expected in 
more massive haloes such as clusters of galaxies (White & 
Frenk 1991; De Lucia, Kauffmann & White 2003). 

We have also tested the impact of a number of assump- 
tions made in our model, and the relevance of numerical 
parameters. Relaxing the IRA for SNII has only a negligi- 
ble effect and does not change the chemical properties of the 
simulated systems. However, as also shown by previous work 
(e.g. Raiteri, Villata & Navarro 1996; Chiappini et al. 1997; 
Mosconi et al. 2001), the chemical properties of both the 
stellar populations and the ISM are sensitive to the time- 
delay associated with SNIa explosions. In agreement with 
previous results, we found that the latter must be included 
in order to be able to reproduce the observed chemical prop- 
erties of galaxies. 

We have shown that chemical enrichment can signif- 
icantly modify the chemical and dynamical properties of 
forming galactic systems. However, the strength of the influ- 
ence of metals on galaxy formation might also depend on the 
physics of energy feedback. If the latter is neglected, as we 
have done here, metals are not mixed widely in gas that has 
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Figure 15. Zero-point of the relation 12+log(0/H) for the gas 
components of the selected galactic objects from the cosmologi- 
cal simulations as a function of the corresponding slopes. We have 
included results from primodial (triangles), metal-dependent (di- 
amonds) and suprasolar (asterisks) cooling functions. For com- 
parison, results for the tests of the isolated disc galaxy run 
with primordial (squared triangle), metal-dependent (squared di- 
amond) and suprasolar (squared asterisk) cooling functions are 
also shown. We have also included observations from Zaritsky et 
al. (1994, filled circles). The diamond and the triangle with the 
arrows have slopes of -0.24 and -0.23, respectively. 



yet to cool. Effects due to metal-dependent cooling remain 
then moderate. In a forthcoming paper, we will complement 
our model for chemical enrichment with a new description 
of a multiphase ISM and energy feedback, which is capa- 
ble of reproducing strong outflows and regulating the star 
formation activity. This combined model should provide a 
powerful tool to exploit observational data on the chemical 
properties of galaxies to constrain galaxy formation. 
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